Mon. Not. R. Astron. Soc. 000,II]-16 (2006) Printed 5 February 2008 (MN WT^ style file v2.2) 



Periodic long-term X-ray and radio variability of 
Cygnus X-1 

Pawel Lachowicz,^* Andrzej A. Zdziarski,^* Alex Schwarzenberg-Czerny/'^* 
Guy G. Pooley^ and Shunji Kitamoto^ 

^Centrum Astronomiczne im. M. Kopernika, Bartycka 18, 00-716 Warszawa, Poland 

Obserwatorium Astronomiczne, Uniwersytet A. Mickiewicza, Sloneczna 36, 60-286 Poznan, Poland 
^ Cavendish Laboratory, J. J. Thomson Avenue, Cambridge CB3 OHE 

'^Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukro, Toshima-ku 171-8520, Japan 



Accepted 2006 July 21. Received 2006 February 14; in original form 2005 August 15 



ABSTRACT 

We present a comprehensive analysis of long-term periodic variability of Cyg X-1 using 
the method of multiharmonic analysis of variance applied to available monitoring data 
since 1969, in X-rays from Vela 5B, Ariel 5, Ginga, CGRO and RXTE satellites 
and in radio from the Ryle and Green Bank telescopes. We confirm a number of 
previously obtained results, and, for the first time, find an orbital modulation at 15 
GHz in the soft state and show the detailed non-sinusoidal shape of that modulation 
in the hard state of both the 15-GHz emission and the X-rays from the RXTE/ PiSM. 
We find the CG'i?0/BATSE data are consistent with the presence of a weak orbital 
modulation, in agreement with its theoretical modelling as due to Compton scattering 
in the companion wind. We then confirm the presence of a ~150-d superorbital period 
in all of the data since ~1976, finding it in particular for the first time in the Ariel 5 
data. Those data sets, covering >65 superorbital cycles, show a remarkable constancy 
of both the period and the phase. On the other hand, we confirm the presence of 
a '^290-d periodicity in the 1969-1979 Vela 5B data, indicating a switch from that 
period to its first harmonic at some time < 1980. We find the superorbital modulation 
is compatible with accretion disc precession. Finally, we find a significant modulation 
in the RXTE/ ASM data at a period of 5.82 d, which corresponds to the beat between 
the orbital and superorbital modulations provided the latter is prograde. 

Key words: accretion, accretion discs - binaries: general - stars: individual: Cyg X-1 
- X-rays: observations - X-rays: stars. 



1 INTRODUCTION 

Cyg X-1, a persistent Galactic black-hole binary, was dis- 
covered in 1964 June (Bowyer et al. 1965), and it has been 
extensively studied since then. Its orbital period is Porb = 
5.6 d and the optical component, HDE 226868/V1357 Cygni 
(Bolton 1972; Webster & Murdin 1972) is an OB supergiant 
(Walborn 1973). The masses of the two stars remain the sub- 
ject of some dispute, with the most recent determination of 
40 ± IOMq and 20 ± 5Mq for the supergiant and the black 
hole, respectively (Ziolkowski 2005). 

Focused wind is the main accretion channel, though 
the companion nearly fills its Roche lobe (Gies & Bolton 
1986; Gies et al. 2003, hereafter G03). Bound- free absorption 
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by the wind leads to a strong modulation with the orbital 
period of the X-rays (Priedhorsky, Brandt & Lund 1995; 
Zhang, Robinson & Cui 1996; Wen et al. 1999, hereafter 
W99; Brocksopp et al. 1999a, hereafter B99; Kitamoto et al. 
2000, hereafter KOO) although the modulation is relatively 
irregular, manifesting itself by dips in the X-ray light curve 
near the superior conjunction of the black hole (Baluciriska- 
Church et al. 2000). The modulation is also seen in the radio 
(Pooley, Fender & Brocksopp 1999, hereafter P99; B99), be- 
ing then due to free-free absorption by the wind (Brocksopp, 
Fender & Pooley 2002). The optical emission is also modu- 
lated with the orbital period, showing two minima/maxima 
per period due to the ellipsoidal shape of the star (Lyutyi 
1985; Kemp 1987; Voloshina, Lyutyi & Tarasov 1997; Brock- 
sopp et al. 1999b). 

Typical of X-ray binaries, Cyg X-1 shows two main 
spectral states, hard and soft (e.g., Gierlihski et al. 1997, 
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1999; Frontera et al. 2001; McConnell et al. 2002; Zdziarski 
& Gierlinski 2004) . Interestingly, the orbital modulation dis- 
appears in the soft state (W99) , probably due to strong pho- 
toionization of the wind by the enhanced soft X-ray flux. 

Cyg X-1 was also found to display superorbital modu- 
lation, i.e., at much longer time scales than Porb- Originally, 
a period of 294 d was reported by Priedhorsky, Terrell & 
Holt (1983, hereafter PTH83) in X-rays and by Kemp et 
al. (1983) in the optical. Later, however, a ~150 d period 
from the radio to X-rays was found by numerous authors 
(B99; P99; KOO; Ozdemir & Demircan 2001; Karitskaya et 
al. 2001; Benlloch et al. 2001, 2004). The probable cause of 
this modulation is precession of the accretion disc. 

We present here a comprehensive analysis of the or- 
bital and superorbital modulation of Cyg X-1 in the X-ray 
and radio bands. We use most of the available X-ray and 
radio monitoring data, applying to it the method of mul- 
tiharmonic analysis of variance (hereafter abbreviated as 
mhAoV). Section|21presents details of our data selection, and 
Section|3 our analysis method. Sections2]and|21present our 
results on the orbital and superorbital modulation, respec- 
tively. Discussion and theoretical interpretation are given in 
Section|Sl and our conclusions, in Section|7| In Appendix A, 
we separately reanalise the presence of the ~290-d periodic- 
ity in the Vela 5B data. In Appendix B, we use the method 
of structure function as an independent check of our results. 



2 DATA SELECTION 

We analyse available light curves separately for the hard 
and soft states. We follow the definitions of those states 
of Zdziarski et al. (2002); thus the state boundaries 
may slightly change if other definitions are used. We 
use X-ray data from the Burst and Transient Source 
Experiment aboard Compton Gamma Ray Observatory 
(CGi?0/BATSE; Paciesas et al. 1997), and the All-Sky 
Monitors (ASM) aboard Rossi X-ray Timing Explorer 
{RXTE; Bradt, Rothschild & Swank 1993; Levine et al. 
1996), Gmga (Makino et al. 1987; Tsunemi et al. 1989), 
Ariel 5 (Holt 1976), and Vela-SB (Conner, Evans & Belian 
1969). We use the 15-GHz radio data from the Ryle Tele- 
scope of the MuUard Radio Astronomy Observatory, and the 
2.25 and 8.30 GHz data from the Green Bank Interferome- 
ter (GBI) of the National Radio Astronomy Observatory in 
Green Bank, WV. Table Ogives the log of the data, which 
are presented graphically in Fig. 

For the RXTE/ KS}A, we use the dwell count rates from 
the Definitive Products database of the High Energy As- 
trophysics Science Archive Research Center (HEASARC). 
Each dwell lasts ~90 s (Levine et al. 1996) and its number 
within one day varies up to 30 for Cyg X-1. For the hard 
state, we have excluded high-flux intervals dominated by the 
so-called failed state transitions (Pottschmidt et al. 2003). 
In this selection procedure, we have chosen MJD 50660- 
50990 as the reference interval, during which the count rate 
in each of the three energy bands (1.5-3, 3-5 and 5-12 keV) 
remained within 4cr around the respective mean. Then, we 
have searched 30-day intervals in the remaining parts of 
the hard state, and have excluded the intervals with more 
than 40 per cent of points in that interval exceeding the 4(t 
level. This, with some further minor adjustments, resulted 



Table 1. The log of the light curves used in this work. H and 
S refer to the hard and soft spectral state, respectively. The two 
sub-intervals of the hard state from the RXTE/ ASM and Ryle 
instruments are treated jointly in the analysis. See Section 1^ for 
the criteria of data selection and for some additionally rejected 
intervals. 



State 


Instrument 


Start 


End 


Time span 






(MJD) 


(MJD) 


(d) 


H 


RXTE/ASM 


50350 


52099 


1750 






52565 


52770 


206 




Ryle 


50377 


52164 


1788 






52565 


52770 


206 




GBI 


50409 


51823 


1415 




CGRO/BATSE 


48371 


51686 


3316 




Ginga/ASM 


46887 


48531 


1645 




Ariel 5/ ASM 


42830 


44292 


1463 




Vela 5B/ASM 


40368 


44042 


3675 


S 


_RXT_E/ASM, Ryle 


52165 


52555 


391 






52800 


52853 


54 



in the intervals of MJD 50590-50660, 50995-51025, 51440- 
51640 and 51840-51960 excluded from the flrst part of the 
hard state of Tabled We have furthermore omitted negative 
count rates, which comprise <0.04 of the ASM data points, 
and thus their omission has a negligible effect. We have also 
applied the barycenter correction to the ASM times of mea- 
surements; however, its effect is negligible. We note that the 
behaviour of Cyg X-1 after MJD 52853 untfl now (~53600) 
has not been suitable for our analysis as showing alternate 
periods of the hard and soft states of rather short durations. 

We notice that there is no agreement in the literature 
on the deflnition of the soft state. For example, Belloni et 
al. (1996) called the 1996 soft state of Cyg X-1 'interme- 
diate' whereas numerous other authors termed it 'soft'. In 
particular, the two occurences of the soft state considered 
by us (see Table Q , while characterized by the relatively 
high and steady high X-ray flux (see Fig.^), may be classi- 
fied as containing also the intermediate state based on some 
other criteria. In general, blackbody-disc dominated states 
of black-hole binaries range between ultrasoft, with almost 
no high-energy tail, to very high or intermediate, where the 
high-energy tail starts at the top of the disc blackbody (see, 
e.g., Zdziarski & Gierlinski 2004 for a review), with differ- 
ent timing properties of the disc blackbody and the tail. The 
soft state of Cyg X-1 lies inbetween the ultrasoft and very 
high states, see, e.g., the spectra in Gierlinski et al. (1999), 
Frontera et al. (2001), Gierlinski & Zdziarski (2003). 

We use the same BATSE data as Zdziarski et al. (2002) . 
They contain 2729 days of usable observations between 1991 
April 25 and 2000 May 22 in the energy bands of 20-100 
keV and 100-300 keV. However, we exclude the periods of 
the soft state of 1994 and 1996, i.e., MJD 49250-49440 and 
MJD 50230-50307. Both of those occurences of the soft state 
are of relatively short duration. For the former, only the 
BATSE data are available, and the latter has already been 
extensively studied, with W99 finding no orbital modulation 
in the RXTE/ASM data. 

We use the Ariel 5/ASM 3-6 keV and Vela SB/ ASM 3- 
12 keV (averaged prior to analysis over 0.25-d intervals) light 
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Figure 1. (a) The RXTiJ/ASM 1-day average light curves corresponding to hard states and the soft states of 2001-2002 and 2003. The 
used boundaries of the soft states are marked by dotted lines. The shaded regions correspond to the data not taken into account in our 
analysis in order to provide approximately the same flux level, (b) The 1-day average 15 GHz light curve from the Ryle telescope, (c) 
The 2.25 and 8.30 GHz GBI data, (d) The CGRO/BATSE light curves. The shaded regions correspond to the 1994 and 1996 soft states, 
which are not included in the analysis, (e) The Ginga/ASM data, (f) The Ariel 5/ASM light curve, (g) The Vela 5B/ASM 0.25-day 
average light curve. 
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curves from the HEASARC database. We note that we use 
the Ariel 5 data (TableEll excluding the MJD 42338-42829, 
which were dominated by soft flare events and included the 
soft state of 1975 (Liang & Nolan 1984). We have also omit- 
ted data points with negative fluxes. We then use the same 
Ginga/ ASM 1-20 keV data as KOO. They contain 339 mea- 
surements during 1987 February-1991 October. 

The Ryle radio data used by us contain 10-min. av- 
erage flux measurements. Parts of this data set have been 
studied by P99 and Benlloch et al. (2004). We have divided 
the data into the hard and soft states (Table based on 
the RXTE/ ASM data. Interestingly, we see a relavitely high 
level of the 15-GHz emission during the soft states. The 2.25 
and 8.30 GHz GBl data (with the average sampling interval 
of ~0.9 d; shown in Fig.Qwith la errors) correspond to the 
hard state only. They have been screened according to the 
standard procedure^ . An early part of the GBI data set has 
been studied by P99 and B99. 



3 TIME ANALYSIS 

In our analysis, we first determine the significance of the 
presence of a period in a light curve using the method of mul- 
tiharmonic analysis of variance (mhAoV; Section 13.111 . We 
use the logarithm of either fiux or count rate, G{t) = \nj^{t), 
similarly to the use of magnitudes in optical astronomy. We 
do not take into account measurement errors in the analysis. 

For a formulation and derivation of the quantitative 
results of the mhAoV method, see Schwarzenberg-Czerny 
(1999) and references therein. Then, for an established pe- 
riod, we fold and average the light curve and fit it with a 
Fourier series, which enables us to find the detailed shape of 
a periodic modulation including its phase fSection l3.21 . 

3.1 The method of period detection 

3.1.1 General formulation and model orthogonality 

In general, we learn from experiments by fitting data, x, 
with a model, a;|| . The data contain n measurements, and the 
model, n|| free parameters. The consistency of the data with 
the model is measured by a function, 0, called a statistic. A 
given model, x\^, using a given statistic (e.g., x'^), yields its 
particular value, Oi. Various methods used in the analysis 
of time series differ both in their choice of the model and the 
statistic; hence are difficult to compare directly. To enable 
such a comparison and for determining the significance of 
results, O is converted into the false alarm probability. Pi. 
This is done considering a hypothetic situation. Hi, in which 
X is pure white noise. Then each pair {x^^,0) corresponds 
to certain cumulative probability distribution of O, namely 
P(n|| , n; O), with Pi being the tail probability that under the 
hypothesis Hi the experiment yields O > 0i, i.e., Pi(0 > 
ei) = l-P(n||,n;ei). 

Up to here, we have just outlined the classical Neyman- 
Pearson procedure of statistics. The specific method for 
analysis of time series used here differs from those commonly 
encountered in astronomy only in the choices of and O. 
Then, our accounting for variance, correlation and multiple 

^ ftp://ftp.gb.nrao.edu/pub/fghigo/gbidata/gdata/OOREADME 



frequencies in calculating P is dictated by the laws of statis- 
tics. The probabilities derived by us from the data are the 
false alarm probabilities. However, we also call them below 
just probabilities or significance levels. 

We note then that Fourier harmonics are not orthogo- 
nal in terms of the scalar product with weights at unevenly 
distributed observations. Certain statistical procedures em- 
ploying classical probability distributions hold for orthog- 
onal models only and fail in other cases. To avoid that, a 
popular variant of the power spectrum, Lomb (1976) and 
Scargle (1982, hereafter LS) periodogram relies on a spe- 
cial choice of phase such that the sine and cosine functions 
become orthogonal. We extend this approach by employing 
Szego orthogonal trigonometric polynomials as model func- 
tions. A series of ny — 2N+1 polynomials corresponds to the 
orthogonal combinations of the A'^ lowest Fourier harmonics 
(Schwarzenberg-Czerny 1996). Orthogonal series are opti- 
mal from the statistical point of view because, by virtue of 
the Fisher lemma (Fisz 1963; Schwarzenberg-Czerny 1998), 
they guarantee the minimum variance of the fit residuals for 
a given model complexity (given by ny). Szego polynomials 
are also convenient in computations since the least-square 
solution may be obtained using recurrence and orthogonal 
projections, resulting in high computational efficiency, with 
the number of steps oc A'^ instead of for N harmonics. 



3.1.2 Variance, the AoV statistics, and model complexity 

The LS method employs the sine as a model, and the 
quadratic norm, 0^2 = ||a; — a;|||p, as the statistic. The cor- 
responding probability distribution is with 2 degrees of 
freedom. Prior to use of the distribution, Q^2 has to 
be divided by the signal variance, V. However, is usually 
not known and has to be estimated from the data them- 
selves. Then, neither 0^2 and variance estimates are inde- 
pendent nor their ratio follows the x distribution, which 
effect has to be accounted for. A simple way to do it 
is to apply the Fisher Analysis of Variance (AoV) statis- 
tic, Q = (n — n||)||a;|| ||^/(n|| ||a; — 2;||||^). Hence we call our 
method, involving Szego polynomials model and the AoV 
statistics, the multi-harmonic analysis of variance or mhAoV 
(Schwarzenberg-Czerny 1996). The probability distribution 
is then the Fisher-Snedecor distribution, F, rather then x^, 
and Pi = 1 — F{n^^,n±;0) where n± = n — ny. For every- 
thing else fixed, replacing x^ with F for n = 100 yields an 
increase of Pi(x^) = 0.001 to Pi(P) = 0.01. Thus, account- 
ing for the unknown variance yields the mhAoV detection 
less significant, but more trustworthy. In this work, n usually 
is larger, for which Pi(P)/Pi(x^) reduces to several. 

Apart from the choice of the statistic, our method for 
TV = 1 differs from the LS one in the average fiux being 
subtracted in the latter (thus yielding ny = 2) whereas a 
constant term is fitted in the former (which can be often of 
significant advantage, see Foster 1995). If the periodic mod- 
ulation in the data differs significantly from a sinusoid (e.g., 
due to dips, eclipses, etc.), then our A'' > 1 models account 
for that more complex shape and perform considerably bet- 
ter then the LS one. For example, we show in Section l4.1l that 
a folded RXTE/ ASM light curve is much better described 
by a model with N — 2 than N = 1 with the probability of 
that improvement being by chance of ~ 10~^^. 
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3.1.3 Multiple trials 

Probability can be assigned to a period found in data ac- 
cording to one of two statistical hypotheses. Namely, (i) 
one knows in advance the trial frequency, vq (from other 
data), and would like to check whether it is also present in 
a given data set or (ii) one searches a whole range, Av, of 
A^efi frequencies and finds the frequency, v, corresponding to 
the most significant modulation. The two cases correspond 
to the probabilities Pi and -Pjv^ff to win in a lottery after 
1 and A^eff trials, respectively, i.e., they represent the false 
alarm probabilities in single and multiple experiments, re- 
spectively. They are related by Pjv„ff = 1 — (1 — Pi)^"^' . Note 
that the hypothesis (ii) and the probability Pjv„ff must be 
always employed in order to claim any new frequency in the 
object under study. The hypothesis (i) is rarely used. How- 
ever, since Pi < Pn^k, it is the more sensitive one. For this 
reason, we advocate its use in the situations where the mod- 
ulation frequency is already known, and we aim at checking 
for its manifestation in the same object but in a new band, 
new data set, etc. We stress that we do not use the hypoth- 
esis (i) to claim any new frequency. 

In this work, we use Pi for cases with the orbital peri- 
odicity and the ~150-d superorbital period, found before in 
numerous studies, and use Pjv^,.f for all other periodicities. 
We calculate Pi and Pjv^ff applying all due corrections (see 
below), and requiring stability of a period after removal of 
a modulation with another frequency (prewhitening) . 

An obstacle hampering use of the (ii) hypothesis is 
that no analytical method is known to calculate A^eff. The 
number A'cff corresponds to independent trials, whereas val- 
ues of periodograms at many frequencies are correlated be- 
cause of the finite width of the peaks. Si/, and because of 
aliasing. As no analytical method is known to determine 
A'eff, Monte Carlo simulations have been used (e.g., Pal- 
tani 2004). Here, we use a simple conservative estimate, 
A'cff = mm{Ai'/5i', Ncaic,n), where A'^caic is the number of 
the values at which the periodogram is calculated. The esti- 
mate is conservative in the sense that it corresponds to the 
upper limit on Pn^h, and thus the minimum significance 
of detection. This effect applies to all methods of period 
search (Horne & Baliunas 1986). In general, it may reduce 
significance of a new frequency detection for large A'cff as 
Pjv^jj S> Pi. In practice, it underscores the role of any prior 
knowledge, in a way similar to the Bayesian statistics: with 
any prior knowledge of the given frequency we are able to 
use the hypothesis (i) to claim the detection with large sig- 
nificance (small Pi). 



3.1.4 Correlation length 

The Pi , and other common probability distributions used to 
set the detection criteria, are derived under the assumption 
of the noise being statistically independent. Often this is not 
the case, as seen, e.g., in light curves of cataclysmic variables 
(CVs). The correlated noise, often termed red noise, obeys 
different probability distribution than the standard Pi , and 
hence may have a profound effect. For example, noise with 
a Gaussian autocorrelation function (ACF) correlated over 
a time interval, St, yields a power spectrum with the Gaus- 
sian shape centered at z/ = and the width = 1/St. It 
may be demonstrated that the net effect of the correlation 



on Pi in analysis of low frequency processes is to decimate 
the number of independent observations by a factor ricorr, 
the average number of observations in the correlation inter- 
val St (Schwarzenberg-Gzerny 1991). Effectively, one should 
use nx/wcorr and S/ricorr instead of nj_ and Q in calculating 
Pi. This result holds generally, for both least squares and 
maximum likelihood analyses of time series. In the respect 
of the red noise, accretion-powered X-ray sources resemble 
CVs and for them the red noise can also have profound con- 
sequences, see e.g., the simulations of Kong et al. (2002). In 
order to take into account this effect here, we check the ACF 
of the fit residuals. 

For independent observations, m — 2 consecutive resid- 
uals have the same sign on average (e.g., Fisz 1963). Thus, 
counting the average length, m, of series of residuals of the 
same sign provides an estimate of the number of consecu- 
tive observations being correlated, ricorr. Note that m — n/l 
where / is the number of such series (both positive and neg- 
ative) . For correlated observations, the average length of se- 
ries with the same sign is m = 2ncorr, which allows us to 
calculate ricorr. 

Let O denote the Fisher-Snedecor statistics from the 
mhAoV periodogram (i.e. from Fourier series fit) computed 
for ny = 2A'-|-1 parameters, n observations and n± = n — ny 
degrees of freedom. To account for ricorr, we calculate Pi as 
follows. 



Pi = 1 - P (^n|| 
where 



e 



tcorr ' tcor 



/_n±_ nil \ 

V2n corr ^ / 



-Hriiie ' 



(1) 



(2) 



and I^ia,})) is the incomplete (regularized) beta func- 
tion (Abramowitz & Stegun 1971), see Schwarzenberg- 
Gzerny 1998 and references therein. (In the popular applica- 
tion, MATHEMATICA, Wolfram 1996, that function is called 
BetaRegularized.) 



3.2 The method of fitting folded Hght curves 

After establishing the existence of a period, we character- 
ize the shape of the modulation by fitting the light curve, 
which we fold and average using a binning appropriate for 
the given statistics and the complexity of its shape. The 
shown errors on the folded and averaged light curves are ob- 
tained from the dispersion (rms) of the fiuxes pre-averaged 
over each real-time bin falling into a given phase bin. (The 
pre-averaging is done because variability over time scales 
shorter than the time length of a phase bin is unimportant 
for the analysis.) In fitting the light curves, we take into ac- 
count the errors obtained in this way. All the uncertainties 
shown below are \a. 

We assume the modulation corresponds to a multiplica- 
tive factor (rather than to an additive fiux component). This 
is suitable to describe modulation due to phase-dependent 
absorption, but it also can describe other effects. Thus, 



T(t) = J^intr(t)Arrod(t), 



(3) 



where Ti^A-cit) is the intrinsic flux and ^mod(i) is the mod- 
ulation factor. 

Our use of the flux logarithm. Git) = \aT(t'), allows 
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us then to describe the modulation as additive. Then we 
use nonlinear least-square method to obtain the best-fitting 
A'^-harmonic Fourier series, 



Gmod(t) = Go - ^ Gfe cos[2Tikv{t - To) - ipk], 



(4) 



where t is measured from some initial moment, To, Go is a 
(fitted) constant term, and we use the same value of A'^ as 
for the corresponding periodogram. The value of A'^ is es- 
tablished by the F-test based on the values of of the fit 
with equation Q. Note that we do not treat v as a, free pa- 
rameter but instead use either the values obtained from the 
periodograms or from other data (e.g., the orbital frequency 
from spectroscopy), compatible with fitting here the folded 
and averaged light curve. 

The obtained modulation can be removed from the light 
curve by, 

Gintr(i) = G(t) — Gmod(t) + Go 



G{t) + ^ Gfc cos[2TTkiy{t - To) 



'Pkl 



(5) 



which preserves the average flux. On the other hand, if the 
modulation is caused by absorption, one can correct the light 
curve for the variable part of the absorption. 



Gintr(t) — G(t) — Gmod{t) + Gn 



(6) 



where Gmax = max[Gmod(i)], which increases the average 
logarithm of the flux by (Gmax — Go), see Tables 2 and 3. 

We remove the orbital modulation prior to searching for 
longer periods (prewhitening) , using equation @, in order 
to reduce the noise level. We then use the periodogram to 
identify the frequency of the next most prominent modu- 
lation, and fit the Fourier series of equation |(1J. Then we 
prewhiten it again, i.e., divide the data by the new modu- 
lation. We repeat the procedure until the peak amplitude 
becomes close to the noise level. To test the persistence and 
significance of the detected modulation, we then repeat the 
procedure using a half of the data. On the other hand, we 
remove the main superorbital modulation found by us (see 
Section 1^ below) before calculating final results for the or- 
bital modulation. 

We define the relative modulation depth by 



^ max ^ min 
max 



(7) 



where J-min and .T-'max are the minimum and maximum fiux 
of the fitted model. This definition, though it differs from 
the standard peak-to-peak one, is suitable for modulation 
caused by absorption. For A'^ = 1, the modulation depth 
and its standard deviation are. 



<^mod = 1 - exp(-2Gi) 



od = 2exp(-2Gi)5Gi, (8) 



where 5Gi is the standard error of Gi. For A'^ > 1, <^niod is 
found numerically, and (5c/>mod, by propagation of errors. 

3.3 Spectral windows 

The discrete spectral window (DSW), 



y ^ cos 2-Kvtk ] + I ^ ^ sin 2-Kvtk, 



(9) 



•aO.02 - 



■^0.01 - 




-6 -5 
log[Frequency (Hz)] 

Figure 2. The spectral window of the RXTE/ ASM light curve. 
The peaks at 2.2 X 10"^ Hz (52.6 d) and 1.2 X 10"^ Hz (0.96 d) 
appear to be caused by the precession of RXTE and observation 
scheduling, respectively. 



where tk is the time of the fcth observation, identifies fea- 
tures related to regularities of the light curve sampling (e.g.. 
Deeming 1975). We have found that only the DSW of the 
RXTE/ ASyi contains interesting features. Namely, there are 
two strong peaks at 52.6 d and 0.96 d and their harmonics, 
see Fig. |5| We have also calculated the DSW for three other 
bright X-ray sources. Crab, Cyg X-3 and GX 5-1, and found 
the same peaks. 

The likely cause of the 52.6-d peak is precession of 
the RXTE satellite, with period in the range of 50-60 d 
(A. M. Levine, personal communication). The precession 
causes a modulation in the frequency of passing through the 
South Atlantic Anomaly, when the detectors are switched 
off. The ~l-d feature appears to be due to RXTE schedul- 
ing (Zdziarski et al. 2004; Farrell, O'NeiU & Sood 2005). A 
peak of the DSW at a frequency v may give rise to aliases of 
actual modulation at frequency offsets of ±i^. However, we 
have found no such features in the studied periodograms. 



4 ORBITAL MODULATION 
4.1 The hard state 

Whenever we find a period of ~ 5.6 d, it is compatible within 
the measurement errors (Table 2) with the spectroscopic 
orbital period. Given this lack of evidence of any difference 
between the period of the modulation and the actual orbital 
period, we use the best available ephemeris for all the data 
sets in fitting the modulation shape. We use the period, V, 
of Brocksopp et al. (1999b) and choose To based on LaSala 
et al. (1998), which gives the ephemeris of 

min[MJD] = 50234.79(±0.01) + 5.599829(±0.00002)£;, (10) 

where E is an integer and To was chosen to be as close as 
possible to all of the intervals of our main data sets. 

Following the method of Section |21 we have found that 
the orbital modulation has a complex shape for some data 
sets, for which the sinusoidal description (A'' = 1) is clearly 
insufficient. Using the F-test (e.g., Bevington & Robinson 
1992) for the fits on the rebinned folded light curves, we have 
found that A' = 3 harmonics of the Fourier series of equation 
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Figure 3. Comparison of fits to the orbital modulation in the 
1.5-3 keV range of the RXTE/ ASM with one (dots) and three 
(solid curve) Fourier harmonics. We see that the single sinusoid 
fails to describe the shape of the modulation. 

(0J are required, in particular, for the RXTE/ KS}A data. 
In the case of the 1.5-3 keV band and for the folded and 
averaged hght curve rebinned into 100 bins, the probability 
that the chance improvement of the fit from = 1 to A*' = 
2 is ~ 10"", and from = 2 to iV = 3, 0.004. Then, 
A'' = 4 is not required, with the probability of 0.67. Fig. |21 
compares the results for A'^ = 1 and A'^ = 3. The shape of 
the orbital modulation can be described as relatively flat for 
the normahzed phase ~0.3-0.7 and with a sharp dip around 
the phase 0/1. 

Our results are given in Table 2 (but note that we 
give the values of n only in Table 3) and Fig. |1] We 
treat jointly the hard-state sub-intervals of the RXTE/ KS}A 
and Ryle data (Table 0. The fractional modulation in the 
RXTE/ ASyi data decreases with the increasing energy, and 
is in good agreement with the results of W99. 

The 15 GHz radio data, see Table 2 and Fig.^fd), also 
require A' = 3. Notably, the minimum of the radio modula- 
tion is significantly shifted with respect to the zero phase, 
by 0.153, or 0.86 d. The corresponding ephemeris (using the 
spectroscopic period) is, 

min[MJD] = 50235. 65(±0. 01) + 5.599829S. (11) 

A similar, but smaller (0.11 or 0.67 d) offset was earlier 
reported by P99. The difference appears to be due to both 
the fit with A' = 1 and a shorter light curve (2 yr) used by 
them. In particular, the fit with A'' = 1 for the current data 
set yields an offset of 0.14. 

The fractional modulation depth of the radio emission 
decreases with the decreasing frequency, see Table 2 and 
Fig.EIe-f). The quality of the 8.3 and 2.25 GHz data is also 
lower than that of the 15 GHz data, and A'' = 1 is sufficient 
to describe the modulation. Interestingly, the offset of the 
minimum increases with the decreasing frequency, reaching 

0. 32±0.09 (1.8±0.5 d) at 2.25 GHz. Since A^ = 1, the values 
of the offset at 8.3 and 2.25 GHz follow directly from Table 

1. These results are similar to those of P99, who studied a 
part of those data up to 1998. 

For the 20-100-keV BATSE data, we calculate the 
mhAoV periodogram with A' = 1, which yields O = 3.7, 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 



Orbital Phase / 2-n 

Figure 4. The phase diagrams of the orbital modulation for the 
RXTE/ ASyi and radio data in the hard state. The solid curves 
represent the fits with N = 3. See Table 2 for details. 



corresponding to the probability of Pi ~ 0.18, see Table 2. 
Although this is a very weak significance, the 20-100 keV 
modulation is probably real given our a priori knowledge 
of the orbital period and the presence of very significant 
modulation in softer X-rays. Our results agree with those 
of Paciesas et al. (1997), who found orbital modulation in 
the 20-300 keV BATSE data. The case for orbital modu- 
lation in the 100-300 keV band is still weaker. With the 
mhAoV method, we still find the period with a small error, 
P = 5.597 ± 0.005, in spite of 9 = 1.0. OveraU, the BATSE 
data are fully compatible with the presence of (weak) mod- 
ulation, but they do not prove it. On the other hand, an 
orbital modulation with parameters similar to those of our 
best fits is implied by the presence of Compton scattering 
in the companion wind, see Section [6.1 1 

The Gtnga/ ASM data show the orbital modulation at 
high significance (Table 2). Our results agree well with those 
of KOO. We note here that the Ginga data clearly show 
higher depths of the modulation than the RXTE/ ASM ones, 
especially for the 6-20 keV band. Taken at the face value, 
the data indicate that the wind from the companion was sig- 
nificantly stronger during 1987-1991 than that after 1996. 

Surprisingly, we find no orbital modulation in the pe- 
riodograms from Vela 5B/ ASM and Ariel 5. This, however, 
agrees with the resuhs of Holt et al. (1979) and PTH83, who 
reported that the 5.6 d modulation was not significantly de- 
tected in the Fourier transforms of those data. 
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Figure 5. The mhAoV periodograms for the 2001-2002 soft 
state showing the lack of evidence for orbital modulation in the 
RXTE/ASM data. 

4.2 The soft state 

We have analyzed two occurrences of the soft state in the 
RXTE/ ASM and Ryle data, 2001-2002 and 2003, see Fig.Cl 
The long duration of the former provides the best available 
constraints. The mhAoV periodograms for the TEXT'S/ ASM 
data are shown in Fig. |3 We see no evidence for orbital 
modulation, confirming the result of W99 for the 1996 soft 
state. This is probably due to a strong increase of the degree 
of wind ionization by the enhanced soft X-ray flux (W99). 

On the other hand, we find orbital modulation in the 
15-GHz data in the soft state of 2001-2002, with V = 
5.6198 ± 0.0043 d at Pi = 0.012 and (/)„iod = (24.5 ± 8.5)%. 
This modulation depth is compatible within errors with that 
in the hard state (Table 2). However, its possible decrease 
can be due to the decrease of the mass loss rate from the sec- 
ondary in the soft state by ~20 per cent, G03 (which also 
contributes to the reduction of the X-ray absoption). The 
relatively low significance level appears to be due to the 
weakness and strong aperiodic variability of the soft-state 
15-GHz flux (Fig.^). Physically, we do expect the modula- 
tion to be present as its mechanism is (most likely) free-free 
absorption (Brocksopp et al. 2002), which is insensitive to 
the ionization level of elements heavier than He (which are 
in turn responsible for the X-ray absorption). 



5 SUPERORBITAL MODULATION 

Here, we search for variability with periods longer than the 
orbital one. We use the same method as in Section 3] How- 
ever, we find that sufficient fits to the folded and averaged 
hght curves can be obtained with N = 1. The main results 
are given in Table 3 and Figs. |S] and |7] 

We have found a ~150 d period to be the most pro- 
nounced superorbital one in the hard-state hght curves. Note 
that we remove the orbital modulation, as described in Sec- 
tion EHK except for Vela 5B/ASM and Ariel 5/ ASM, which 
do not show orbital modulation), and other superorbital 
modulations (whenever present, see the bottom part of Ta- 
ble 3) prior ot calculating the final results for the ~150 d 
superorbital modulation. 



293" 192" 151" 




0.005 0.0075 0.01 0.0125 0.015 



Frequency (d"') 

Figure 6. The low-frequency hard-state mhAoV periodograms 
(after removing the orbital modulation), (a) The RXTE/ ASM 
data, (b) The aCRO/BATSE data, (c) The Ryle radio data, (d) 
The GBI radio data. 

We then determine the common ephemeris to all the 
data using the weighted average, 151.43 ± 0.20 d, of all the 
~150 d periods, and choose To (around MJD 50500) based 
on the RXTE/ ASM 5-12 keV data, which provide the high- 
est signal-to- noise ratio of the modulation (see Table 3). This 
yields 

min[MJD] = 50514.59 + 151.43(±0.20)-E. (12) 

We find that each of the individual periods is consistent with 
the weighted average within < 2g. 

Apart from the V ~150 d, we find also some longer 
periods, >200 d, in the earliest observations, i.e., by Vela 5B, 
Ariel 5, Ginga. For those, we have arbitrarily chosen To = 
44450.0 (MJD). We separately discuss the ~290-d period 
from the Vela 55/ ASM 10-yr monitoring in Appendix 1X1 

Fig. inia) shows the hard-state RXTE/ ASM peri- 
odograms. The most significant peak in all three energy 
bands corresponds to V — 150 d. Similar values have been 
earlier reported in the 1.5-12 keV ASM data by B99, Ben- 
Uoch et al. (2001, 2004) and Ozdemir & Demircan (2001) 
and in optical data by Karitskaya et al. (2001). Also, that 
periodicity is confirmed by the structure-function analysis, 
see Appendix B. We find the fractional modulation in the 
three bands to be compatible with constant, see Table 3. 

Prewhitening with the ~150-d modulation has re- 
sulted in the lack of other significant periodicities in the 
RXTE/ ASM data. Although we see two peaks with longer V 
in the inhAoV periodograms, Fig.|SJa), their position change 
after prewhitening with the orbital and the ~150-d periods. 
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Figure 7. The phase diagrams of the superorbital modulation in the hard state for the ephemeris of equation 1121 . modelled by sinusoids. 
Note the approximate constancy of the phase for all the light curves. See Table 3 for the parameters. 



Also, those periods change when we consider shorter inter- 
vals of the ASM data. In particular, we do not confirm here 
the 182. 5-d period claimed by Ozdemir & Demircan (2001). 

Fig. Etb) shows the 20-100 keV and 100-300 keV 
mhAoV periodograms corresponding to the hard state ob- 
served by BATSE. The only significant peak corresponds to 
V ~ 150 d. It is also confirmed by the structure-function 
analysis, see Appendix B. Our findings contrast those of 
Paciesas et al. (1997), who reported no periodicities at fre- 
quencies ^0.1 d^^. Still, one can see a Fourier peak corre- 
sponding to ^ 150 d in their power spectrum. Then, B99 
found V ~ 142.0 ± 7.1 d in the 20-100 keV data before the 
1996 state transition, in agreement with our findings. Inter- 
estingly, the depth of the modulation is lower than that of 
the 1.5-12 keV emission, see Table 3. 

The 15-GHz hard-state data show V = 150.7 ± 0.6 d, 
see Fig. |SJc). However, that peak is accompanied by an- 
other statistically significant periodicity at 192.3 ± 1.0 d. 



We have checked that the ~ 150 d periodicity remains after 
prewhitening with the V ~ 192 d period. The 15-GHz data 
also show V — 293 d period, similar to the ~294 d X-ray 
modulation found by PTH83 in the Vela 5B data. However, 
prewhitening of the radio hght curve with both the ~150 
d and and 192 d periods causes the change of that longest 
period from 293 d to ~350 d, which casts doubt on the re- 
ality of that modulation. In fact, neither the 192 nor 293-d 
periodicity appear applying the structure-function analysis, 
see Appendix B. 

The GBI 8.3 and 2.25 GHz data also show prominent 
modulation at ~150 d, see Table 3 and Figs. I^Jd) and |7| 
Their depth is compatible with being equal to that of the 
15 GHz emission. Interestingly, the periodograms also show 
peaks around ~190 d, similar to those seen in the Ryle data. 
Subsets of the Ryle and GBI data up to 1998 were earlier 
analyzed by B99, who obtained results compatible with ours. 

The Gtnga/ ASM data were earlier analyzed by KOO. 
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They reported two periods, 150 d and ~ 210-230 d. Here, 
we confirm their results, obtaining two peaks at 151 d and 
230 d. 

In the Ariel 5 3-6 keV hard-state data, we have found 
two prominent peaks (not shown here). The first one &tV — 
278.6 ±3.1 d is also confirmed by the structure-function 
analysis (Appendix B). It is close to the V ~ 296 ± 10 d 
reported by PTH83. 

The second peak corresponds to V — 151.3 ± 0.8 d, 
which remains significant after prewhitening with the first 
period. We note that the ~ 150 d periodicity is seen, in fact, 
in the Fourier spectrum of PTH83. 

Remarkably, all the ~150-d modulations since 1976 are 
consistent with not only having the same period, but also 
show an almost constant phase, see Table 3 and Fig. |7| The 
phase constancy also support our determination of the pe- 
riod from the weighted average, as an error on that would 
result in phase shifts in data separated by long intervals. 

The approximately constant phase is kept in spite of a 
number of soft-state occurrences interrupting the hard state. 
The longest available soft-state period, of 2001-2002 is still 
too short (391 d. Table to enable us to study the super- 
orbital modulation in that state itself. 



6 DISCUSSION 

6.1 Orbital modulation 

We consider the usual spherically symmetric approximation 
to the radial density profile of the wind from the companion. 



n{R) 



no 



[1 - {R./R)Y 



(13) 



(e.g.. Castor, Abbott , 
the velocity profile of 

R 



v{R) 



Klein 1975), which corresponds to 



(14) 



and where i?* is 

^2 



the radius of the star, no ~ 
Mioss/(rnii4nRiva:,) is the density parameter, Mioss is the 
mass-loss rate, mn is the hydrogen mass, Vaa is the terminal 
velocity of the wind, and a is the power-law index of the 
dependence. 

These formulae were used by W99, who have success- 
fully explained the orbital modulation seen in the hard state 
in their RXTE/ ASM data by phase-dependent absorption in 
the ionized wind. They found i ~ 30^20° as the inclination 
implied by the modulation profile. 

Here, we consider modulation due to Thomson scatter- 
ing by the wind, relevant to the BATSE 20-100 keV energy 
range. We obtain then the fractional modulation depth [eq. 
GJ] of 



exp[-rT(7r)] - exp[-rT(0)] 
exp[— rT(7r)] 



(15) 



where tt is the Thomson optical depth integrated along the 
line of sight from the black hole through the wind. 



(16) 



TTlt/Jorbj = CTT / n[R{r,(t>oih,i)]dr 

'0 

and ctt is Thomson cross section. We have assumed the 
rate of mass loss by the companion during the hard state 




30 

Inclination 

Figure 8. Tiie fractional orbital modulation due to Thomson 
scattering in the wind as the function of the inclination, see Sec- 
tion |0] 



of Mioss 1.6 X 102°gs-i (G03), 7?* ~ 1.58 x 10^^ cm, 
the binary separation, a = {P /2-n-)^^^[G{Mi + M2)Y^^, of 
2.27i?* (Ziolkowski 2005), and Woo ~ 1586kms~\ a = 1.05 
(Gies & Bolton 1986). Then, no ~ 1.9 x 10^" cm'^ The 
results are shown in Fig. |H| The observed orbital modula- 
tion in the Thomson regime (20-100 keV) of ~3 per cent 
corresponds to i ^ 35° , and is in agreement with the corre- 
sponding result of W99. These theoretical results, based on 
optical data, support the reality of the orbital modulation 
in the BATSE data. We caution, however, that the actual 
wind in Cyg X-1 is clearly not spherically symmetric, and 
thus those determinations of i bear an additional systematic 



The radio modulation is due to free-free absorption of 
the jet emission in the wind (Brocksopp et al. 2002^). The 
phase shift of the modulation with respect to the orbital 
phase is likely due to the time it takes for the jet material, 
ejected close to the black hole, to propagate to the place 
where the radio emission originates (i.e., jet bending). 



6.2 Superorbital modulation 

A generally accepted interpretation for the superorbital 
modulation, observed also in a number of other binary X- 
ray sources (e.g., SS 433, Her X-1, LMC X-3, LMC X-4), is 
accretion disc precession (see e.g., Wijers & Pringle 1999 for 
a review and references). 

If, due to some reason, the disc is inclined with respect 
to the orbital plane (see Fig. |5J , it will retrogradely precess 
due to the tidal forces exerted by the secondary (Katz 1973). 
This appears to explain the superorbital periodicities in Her 
X-1, SS 433 and LMC X-4 (Gerend & Boynton 1976; Lei- 
bowitz 1984; Heemskerk & van Paradijs 1989). The period 
of the tidally-forced precession is given by. 



3 ^PRl 
7\ a 



3/2 



■ cos 5, 



(17) 



where fi is the ratio of the mass of the secondary to that 
of the compact object, i?L is the Roche lobe radius of the 

Note, however, apparent errors in their equations (2), (6), (8). 
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Figure 9. Schematic view of Cyg X-1 at the zero orbital phase. 
As before, i is the orbital inclination, and & is the angle between 
the inclined disc and the orbital plane. 



primary, /3 is the ratio of the outer disc radius to _Rl, and 

5 is the inchnation of the disc with respect to the orbital 
plane (Larwood 1998). Equation 1171 assumes that the disc 
precesses as a rigid body, which is the case when the sound 
crossing time in the disc is much less than the precession 
time scale (Papaloizou & Terquem 1995). We note that the 
numerical coefficient above depends on the assumed disc 
model, and may be different from 3/7 (e.g., Larwood 1997). 

The Roche lobe radius is well approximated by, 

0-49 .^o^ 

a ~ 0.6 + m''/^ ln(l + At-^/3) ^ ' 

(Eggleton 1983). At 2 < < 3 estimated for Cyg X-1 (Gies 

6 Bohon 1986; G03; Ziolkowski 2005), Ri^/a ~ 0.32-0.29. 
Then, 

1/3 



i?L ^ 10' 



Ml + hh 



50 M, 







(19) 



In the case of accretion via Roche-lobe overflow and for 
thin and weakly viscous discs (Paczyriski 1977), /3 ~ 0.68- 
0.66 at 2 < /i < 3. Then for viscous Roche-lobe flows, the 
size of the disc corresponds to the tidal radius, for which 
/? ~ 0.87 (Papaloizou & Pringle 1977). However, Cyg X- 
1 accretes via focused wind, in which case the disc size is 
usually smaller (e.g., Frank, King & Raine 2002). Therefore, 
equation 11711 puts rather weak constraints on the amplitude 
of precession, 5. At = 2, 5 ~ 0°-60° for /3 ~ 0.55-0.87, 
and at ^ = 3, 5 ~ 0°-63° for (5 ~ 0.52-0.87. StiU, equation 
11711 shows that the precession in Cyg X-1 is fully consistent 
with being driven by the tidal force of the companion. 

On the other hand, Wijers & Pringle (1999) have con- 
sidered precession due to radiation-induced warping in X-ray 
binaries. This can yield either retrograde or prograde preces- 



sion. The corresponding precession period is roughly given 
by (see their eq. 18), 



3/4 



M 



Ml 



WMq J \^10i7gs 



-3/10 



where M is the accretion rate, a is the viscosity parameter 
and e is the accretion efficiency. For Cyg X-1, 

(F) 



M ~ 2 X 10' 
d 



X 



2kpc 



4 X 10~* erg cm-^ g-i 

(ot) sb"' 



(21) 



where (F) is the average bolometric hard-state flux and d is 
the distance, and we used above their values of Zdziarski et 
al. (2002) and Ziolkowski (2005), respectively. 

The above estimates yield Pwarp — 10^ d, much longer 
than that observed. We note that Wijers & Pringle (1999) 
have given an estimate of Pwarp — 180 d for Cyg X-1, but 
assuming a and M as high as 1 and 1.4 x lO'^^gs"^, respec- 
tively. Thus, we consider radiation warping as an unlikely 
cause of the precession of Cyg X-1. Furthermore, the ob- 
served precession has a remarkably stable period, changing 
at most weakly over > 15 years, whereas Wijers & Pringle 
(1999) note that their mechanism is unlikely to give a stable 
period for sources with variable luminosity. 

Regardless of the mechanism causing the precession, 
there is the issue of the mechanism causing the flux modula- 
tion during the precession. In principle, there are a number 
of possibilities. One is that the outer edge of the disc (com- 
pletely optically thick) partially covers the X-ray source. 
This, however, would require extreme fine-tuning to achieve 
a ~25 per cent depth of the modulation. Namely, the X-ray 
source has the size ~ 10^ -Rg (where _Rg = GM/c?), as indi- 
cated by the X-ray power spectrum extending to high fre- 
quencies and agreement with theoretical prediction on the 
range of radii where most of the accretion power is released. 
On the other hand, the outer edge of the disc is at a much 
larger distance, > 10^ -Rg (see the discussion above). Another 
possibility is that the outer part of the disc fully obscures 
the X-ray source, but we see X-rays scattered in a large-size 
corona above the disc. This, however, would dramatically 
affect the X-ray power spectrum, resulting in a cutoff above 
~0.1 Hz, which is clearly not seen. Then, bound-free ab- 
sorption in a spatially-extended, moderately optically-thin, 
medium associated with the outer regions of the disc appears 
to be ruled out as there is rather weak energy dependence 
of the modulation. On the other hand, a viable scenario is 
the outer disc wind/corona being almost fully ionized, with 
scattering away from the line of sight being responsible for 
the superorbital modulation. 

Yet another possibility is that the X-ray emission is in- 
trinsically anisotropic. Such a possibility was considered by 
B99, who relied in calculating the implied inclination on the 
blackbody-type anisotropy of a slab with the ffux propor- 
tional to the projected area. However, there is overwhelm- 
ing evidence that the dominant radiative process produc- 
ing X-rays in the hard state of Cyg X-1 (and other black- 
hole binaries) is thermal Comptonization (e.g., Zdziarski & 
Gierlinski 2004). The anisotropy of thermal Comptonization 
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(Poutanen & Svensson 1996) in planar geometries, though 
substantially weaker than that ol the constant specific in- 
tensity case, appears fully capable to explain the observed 
superorbital modulation (work in preparation). 

We note that our result on the stability of the super- 
orbital period over > 15 yr rules out the model for strong 
outbursts of Cyg X-1 by Romero, Kaufman Bernando & 
Mirabel (2002). Those outbursts were observed (Stern, Be- 
loborodov & Poutanen 2001; Golenetskii et al. 2002, 2003) 
on MJD 49727, 49801, 51287, 51289, 52329, 52364 and 
52682. Romero et al. (2002) have proposed that the events 
are caused by precession of the jet of Cyg X-1, with the 
flares corresponding to the strongly beamed emission of the 
jet crossing our line of sight every precession period. They 
considered, in particular, the flares on MJD 49727, 49801, 
52329, reported by Golenetskii et al. (2002). The first two 
are separated by 74 days, which is about a half of the pre- 
cession period of Cyg X-1. Romero et al. (2002) proposed 
that that precession period undergoes occasional changes, 
which would explain the value of 74 d. However, our present 
results show a remarkable stability of the precession period 
over the last several tens of years and rule out that expla- 
nation. Furthermore, the larger number of the dates of the 
outbursts compiled by Golenetskii et al. (2003) do not show 
any regularity. Interestingly, the event on MJD 52329 took 
place at the minimum of ~150 d superorbital cycle. 

With the RXTE/kSM data, we have also looked for 
signals corresponding to the beat between the orbital and 
superorbital frequencies. If the physical cause of the super- 
orbital frequency is modulation of the X-ray flux emitted in 
a given direction by disc precession, the photons reflected 
(e.g., Magdziarz & Zdziarski 1995) from the surface of the 
secondary will change at the beat frequency of the two mod- 
ulations. Then, the position of a possible peak in the peri- 
odogram at either VoA + i^sup or Vovh — Va^ip would tell us 
whether the precession is retrograde or prograde. We have 
thus searched for the corresponding peaks at the periods of 
5.40 d and 5.82 d. 

Fig. llOf a) shows the periodograms for the 1.5-3 keV 
band before and after prewhitening with both the orbital 
(calculated for TV = 3) and superorbital frequencies. After 
the prewhitening, the use of A'^ = 1 is sufBcient to account 
for the shape of any remaining features. We clearly see a pro- 
nounced peak at 5.82 d. Fig. llOf b') shows that such peaks are 
present in all three ASM bands. The statistical signiflcance 
of the features, Pn^h, and the relative modulation depth, 
4>mod, for the 1.5-3, 3-5 and 5-12 keV bands is 1 x 10~^ and 
0.05, 2.2 X 10"* and 0.03, 2 x 10"^ and 0.01, respectively. 
No statistically significant 5.40-d variability is detected. On 
the other hand, there is a rather strong peak (though much 
weaker than that at 5.82 d) at ~5.5 d in the 1.5-3 keV band, 
which origin is unclear. 

This may indicate that the precession in Cyg X-1 is 
prograde. However, the interpretation of the beat frequency 
as due to variable reflection from the companion surface 
presents some problems. Taking the determination of the 
binary parameters of Ziolkowski (2005), we find that the 
companion subtends a solid angle of ~ 0.05 x in as seen 
from the black hole. The observed modulation amplitude of 
~0.05 in the 1.5-3 keV band would require a unit albedo and 
such precession amplitude that the X-ray source seen from 
the star surface is fully obscured by the disc. Furthermore, 
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Figure 10. (a) The mhAoV periodograms for the RXTE/ ASM 
1.5—3 keV data after prewhitening with both the orbital and su- 
perorbital periods (blue curve) and before it (black curve). Both 
periodograms show peaks corresponding to the period of 5.82 d, 
possibly indicating the prograde character of the precession, (b) 
Detailed shape of the periodograms after the prewhitening in the 
three ASM bands. 



the albedo has to decrease by a factor of several in the 5-12 
keV band. We have also checked the relative phase between 
the orbital and beat modulations. In the above interpreta- 
tion, the two modulations should be in phase around the 
minimum of the superorbital modulation, when the disc is 
most inclined towards the companion. However, we find the 
ephemeris of the beat modulation to be. 



min[MJD] = 50239.70 -I- 5.82E, 



(22) 



which corresponds to the orbital and beat modulations being 
in phase instead around the midpoint between the minimum 
and the maximum of the superorbital modulation. Such rel- 
ative phases appear to have no geometrical interpretation. 

We have also searched for signatures of nutation. This 
effect is due to a perturbation of the disc rotation at the 
maximum of the torque exerted on the disc by the compan- 
ion. This takes place when the star is closest to the point on 
the outer edge of the disc at the maximum distance from the 
orbital plane. Therefore, the nutation frequency is twice the 
beat frequency. We have not found signatures of nutation in 
any of our data. 

We would like to stress the remarkable stability of 
the ~150-d modulation. This modulation is clearly present 
starting with the Ariel 5 data used by us, i.e., the begin- 
ning of 1976, and remains until at least mid 2003. This cor- 
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responds to about 65 precession periods. Furthermore, the 
phase of the modulation has also remained approximately 
constant, see Table 3. This appears to present a problem for 
the interpretation of the prograde character of precession. 
Precession can be prograde due to radiation-induced warp- 
ing (Wijers & Pringle 1999), but it is then unlikely to have 
such constant period given the variable intrinsic luminosity 
of Cyg X-1 and its state transitions. Superhump-type pre- 
cession is prograde, but it can take place only for /i <g 1 
(e.g., Frank et al. 2002). Then, our results may indicate the 
presence of a new mechanism causing stable prograde pre- 
cession at /X > 1. 

On the other hand, the ~150-d period is not present 
in the Vela 5b data, ~1969-1979, where, instead, the ~290 
d period appears (Appendix^ PTH83). Interestingly, we 
may see a slow evolution of this period in a following time 
interval, with secondary superorbital periods in the Ariel 5 
and Ginga data of 278 d and 230 d, respectively. (The change 
of the long period from ~290 d to a half of this value was 
noticed by B99.) 



7 CONCLUSIONS 

We have studied periodic long-term variability of Cyg X-1 
using a statistical method different from those used before, 
and apply it, for the first time, to the data set spanning over 
~30 years. In our method, we calculate the periodograms 
and probabilities with the multi-harmonic analysis of vari- 
ance, and then fit the folded and averaged light curves. We 
have analyzed X-ray monitoring data from Vela 5B, Ariel 5, 
Ginga, CGRO and RXTE, and radio data from the Ryle and 
Green Bank telescopes. We have confirmed and refined pre- 
vious results on the 5.6-d orbital modulation of X-ray and 
radio emission in the hard state, caused by the attenuation 
in the wind of the companion supergiant being dependent 
the binary orbital phase. In particular, we show the detailed 
non-sinusoidal shape of the orbital modulation of X-rays ob- 
served by the RXTE/ ASM and of the 15 GHz emission. For 
the first time, we find an orbital modulation at 15 GHz in 
the soft state. We confirm the presence of orbital modulation 
in the 20-100 keV GGRO/BATSE data. 

Then, we find the presence of a ~150-d superorbital pe- 
riod in all of the data since ~1976; in particular, we find its 
presence for the first time in the Ariel 5 data. Very remark- 
ably, we find both that period and its phase have been stable 
over > 65 superorbital cycles. This coherence poses severe 
restrictions on any underlying clock mechanism. We find the 
cause of the superorbital modulation to be compatible with 
accretion disc precession, probably due to the tidal forces 
exerted by the secondary. Furthermore, we find modulation 
at the frequency corresponding to the beat between the or- 
bital modulation and the prograde disc precession. On the 
other hand, we confirm the presence of a ~290-d periodicity 
in the 1969-1979 Vela 5B data. This indicates a switch from 
that precession period to its first harmonic around ~1980. 
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APPENDIX A: THE ~290-DAY PERIODICITY 
IN THE VELA 55/ ASM DATA 

The earliest reported long-term periodicity found in Cyg X- 
1 data is the 294 ± 4 d one found by Vela 5S/ASM in the 
3-12 keV range (PTH83). The same period was shortly after 
found in optical photometric data by Kemp et al. (1983). In 
Sectional we have found a similar period of 279 ± 3 d in the 
3-6- keV Ariel 5/ ASM data, as well as a tentative ~ 293-d 
period in the 15 GHz radio data. However, that period has 
not been confirmed in any other X-ray data. 

Here, we analyze the Vela 55/ ASM observations of Cyg 
X-1 binned to 0.25 day averages using the mhAoV method. 
The resulting periodogram is shown in Fig. lAll We find a 
very strong dominant peak at the period of 289. 5± 1.7 d and 
the absence of any significant feature at the first harmonic. 
The significance, phase, and the depth of the modulation 
are given in Table 3. We have also rebinned the data to 30-d 
bins, as done by PTH83, and found the maximum slightly 
shifts to 294 ± 3 d. Thus, our results are in complete agree- 
ment with those of PTH83. 

We should note that the precession period of Vela 5B 
is --300 d, see the Vela 5B/ASM Calibration Guide of 
HEASARC ^. This could, in principle, suggest an instrumen- 
tal origin of the observed periodicity. In order to test this 
possibility, we have obtained Vela 5B/ ASM periodograms 
for three other bright X-ray sources, 4U 0614+09, Cyg X-3 
and Her X-1, using the same time span (see Table 1), energy 
band and binning as for Cyg X-1. In all of those sources, we 
have found no significant periodicities around ~300 d. These 
results as well as the confirmation of the ~290-d period in 
the optical data (Kemp et al. 1983) appear to confirm the 
reality of the peak in the Vela 5B/ASM periodogram. 

To further test the origin of the observed periodicity, 
we have generated a Monte- Carlo light curve assuming the 

heasarc. gsfc. nasa.gov / docs/heasarc / caldb /docs/ vela5b 
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count rate equal to the average one of Cyg X-1 with a Gaus- 
sian distribution of the errors at the signal-to-noise ratio 
matching the data at the times matching those of the obser- 
vations. We have then calculated the mhAoV periodogram 
and found no significant ~300-d or other flux periodicity. 
In fact, < 5 was found in the 0.003-0.2 d~^ frequency 
range. This further confirms the reality of the presence of 
the ~290-d period in the data. 

We have also examined the periodograms for the time 
interval when Cyg X-1 was observed simultaneously by Arid 
5/ ASM and Vela 55/ ASM, i.e., MJD 42830-44042. We have 
found those periodograms have the strongest peaks at the 
periods of 302 d and 291 d, respectively. There is a ~150-d 
peak in the Ariel 5 data, but only as the third strongest 
one, and no corresponding peak is present in the Vela 5B 
data. This further supports our findings above. 

It then appears that the superorbital period of Cyg X-1 
had indeed changed over a time scale > 10 yr from ~290-d 
to approximately its first harmonic, ~150 d, seen in all the 
observations after those by Vela 5B. During a transitionary 
time interval, both periods were present, see Table 3. 



APPENDIX B: 
OF CYG X-1 



THE STRUCTURE FUNCTION 



A method to quantify long-term variability alternative to 
periodograms is to calculate the structure function (SF). It 
was introduced by Kolomogorov (1941a, b), and applied for 
the first time in astronomy by Simonetti, Cordes & Heeschen 
(1985), and then used by many authors, e.g., Hughes, AUer 
& AUer (1992), Collier & Peterson (2001), and Czerny et al. 
(2003). 

It is a function of the time domain, with the first-order 
SF defined by 



SF{r) = {[F(t + r)-F{t)f), 



(Bl) 



where F{t) is the light curve. If F(t) = A smt, i.e., it is 
sinusoidal with the 2n period and we measure it between 
the times a and h, the structure function can be calculated 
to be. 



SF(r) = 



b — T — a 

2 T 



2 sin 



1 



[sin(t + t) — sint]^ At 
2 cos(a -I- b) sin(a — 6 + r) 



(B2) 
(B3) 



where r < 6 — a, and the second (boundary) term in brackets 
is negligible for r <C fe — a. We see that the structure function 
has minima (at 0) for r equal to the period and its subhar- 
monics, i.e., t = 2fc7r, where k is an integer. This should be 
bear in mind when analyzing results of application of this 
technique to observed light curves. 

We then search for long-term periodicities in our light 
curves corresponding to the hard state (except that of Gtnga, 
due to its limited quality). We use the procedure described 
in Czerny et al. (2003) for calculating the SF and its uncer- 
tainty. The results are shown in Fig. IBll 

We clearly see the dips to the ~150-d period and its 
subharmonics in the RXTE/ ASM, BATSE and Ryle data. 
In the ASM data, the dip position is at 158±14 d. (Hereafter, 
the uncertainty has been estimated as equal to the FWHM 



7.5 

£ 5 
t. 

" 2.5 



40 

30 

h 

fcT 20 
10 


0.05 
S 0.04 
" 0.03 

0.02 

275 
„ 250 
i 225 
m 200 
175 
150 




1.5-3 keV RXTE 
3-5 keV 
5-12 keV 

, _ , H ,in."i"".' 



-+- 



H h- 1- 



100-300 keV 



H 1 1 I I I I I I 

, .*.>.;j;4f': 



Ryle 15 GHz 



3-6 keV Ariel 5 



3-12 keV Vela 5B 



1' S J 



10 100 1000 

T(d) 

Figure Bl. The structure function for our hard-state light 
curves. Periodic variability in the data manifests itself as dips 
at T equal to the period and its subharmonics. 



of the absolute value of the dip profile after subtracting the 
surrounding continuum.) Thus, this result is consistent with 
our periodogram results. In the 20-100 keV BATSE data, 
the dip is at 153 ± 16 d, while the 100-300 keV band shows 
only a shallow minimum at ~150 d. The radio data show the 
dip at 151 ± 19 d. Interestingly, there is no sign of any dip at 
~192 d, which was found in the corresponding periodogram 
(see Section |SJ. On the other hand, the Ariel 5/ ASM show 
a pronounced dip at ~276 d. The Vela SB show only a very 
shallow minimum around ~300 d. 

Thus, the SF technique fully confirms the periodogram 
results for the RXTE/ ASM and for the 20-100 keV BATSE 
band. Thus, we can consider those results as beyond any 
reasonable doubt. For other data, the SF results are in only 
partial agreement with the periodogram ones. 

This paper has been typeset from a TJ^X/ ETJ5X file prepared 
by the author. 
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Table 2. The results of the analysis of the orbital modulation in the hard state. See Section 3 for the definitions. The values of P, 0, ricorr and Pi are found directly from the 
pcriodograms whereas the remaining quantities are obtained for the folded and averaged lightcurves using the ephemeris of equation (10). See Table 3 for the values of n. The units of 
exp(Go) are those shown in Fig. 1. 
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-0.089±0.235 






0.016 


3.1±0 


CGflO(100-300kcV) 


5.5974±0.0049 


1.0 


3.05 


0.81 


1.836±0.004 


0.010±0.008 






-0.266±0.506 






0.010 


1.6±1 


Ginga(l-6kcV) 


5.6014±0.0023 


12.6 


1.24 


2 X 10-*^ 


-0.688±0.016 


0.127±0.024 






0.039±0.173 






0.127 


22.4±3 


Ginga(6-20 kcV) 


5.6027±0.0032 


4.9 


1.09 


4 X 10-^ 


-1.148±0.018 


0.086±0.027 






0.246±0.292 






0.086 


15.8±4 



Table 3. The most significant superorbital periodicities found in the hard state. Af = 1 in all cases, and n gives the total number of points in the light curves. See Appendix A for 
discussion of the Vela 5B data. The values of V, 0, ncorr, Pi and fjv^fj are found directly from the periodograms whereas the remaining quantities are obtained for the folded and 
averaged lightcurves with the ephemeris of equation (12) used for the first 11 rows. Since A'^ = 1 throughout, Gmax — Go = Gi. The units of exp(Go) are those shown in Fig. 1. 



Detector 


r [d] 


© 


n 


TT-corr 


Pi 




A^cff 




Go 


Gi 




</'mod[%] 


RXTE{l.5-3 keV) 


152.60±1.20 


722.7 


27232 


1.80 


< 10- 


10 






1.771±0.015i 


0.134±0.02li 


0.287±0.169i 


23.5±3.2 


RXTE{3-5 keV) 


153.02±0.79 


1233.9 


27344 


1.89 


< 10- 


10 






1.800±0.012i 


0.115±0.016i 


0.124±0.148i 


20.5±2.5 


RXTE{5-12 keV) 


153.10±0.66 


1563.6 


27345 


1.91 


< 10- 


10 






2.182±0.009i 


0.114±0.013i 


0.000±0.122i 


20.4±2.1 


Ryle(15 GHz) 


150.69±0.55 


257.30 


12493 


2.14 


< 10- 


10 






-4.454±0.008i 


0.106±0.012i 


-0.041±0.119i 


19.1±1.9 


GBI(8.30 GHz) 


147.30±1.03 


44.1 


1547 


1.12 


< 10- 


10 






-4.184±0.008i 


0.124±0.013i 


0.325±0.09li 


22.0±2.0 


GBI(2.25 GHz) 


150.92±0.59 


37.3 


1547 


1.18 


< 10- 


10 






-4.247±0.01li 


0.107±0.014i 


0.329±0.160i 


19.3±2.2 


CGi?O(20-100 keV) 


151.04±0.32 


109.1 


2511 


2.28 


< 10- 


10 






2.324±0.005i 


0.081±0.007i 


-0.003±0.089i 


14.9±1.2 


CGi?O(100-300 keV) 


152.37±0.51 


41.9 


2511 


3.05 


< 10- 


10 






1.843±0.007i 


0.070±0.008i 


0.096±0.14li 


13.1±1.4 


Ginga{l-6 keV) 


152.13±1.55 


11.3 


337 


1.22 


7 X 10 


-6 






-0.717±0.017i 


0.109±0.024i 


-0.005±0.223i 


19.6±3.9 


Ginga{6-20 keV) 


151.01±1.52 


10.1 


337 


1.19 


2 X 10 


-5 






-1.179±0.015i 


0.109±0.024i 


0.157±0.208i 


19.6±3.9 


Ariel 5(3-6 keV) 


151.33±0.79 


24.8 


1972 


1.30 


< 10- 


10 






-1.053±0.01li 


0.074±0.015i 


0.477±0.220i 


13.8±2.6 


Ryle(15 GHz) 


192.29±1.04 


281.8 


12493 


2.12 


< 10- 


10 


33.1 


< 10-10 


-4.461±0.01l2 


0.107±0.01l2 


-0.406±0.10l2 


19.3±1.8 


GBI(2.25 GHz) 


185.47±1.10 


27.9 


1547 


1.19 


< 10- 


10 


17.8 


< 10-10 


-4.257±0.0172 


0.052±0.0242 


0.559±0.4462 


9.9±4.3 


Ginga{l-6 keV) 


230.30±3.39 


17.3 


337 


1.26 


2 X 10 


-8 


21.4 


4 X 10-^ 


-0.701±0.0183 


0.106±0.0263 


1.191±0.2363 


19.1±4.2 


Ginga{6-20 keV) 


232.65±3.35 


11.8 


337 


1.11 


1 X 10 


-6 


20.9 


2 X 10-5 


-1.139±0.0143 


0.042±0.0193 


0.615±0.5173 


8.1±3.5 


Ariel 5(3-6 keV) 


278.55±3.07 


81.7 


1972 


1.31 


< 10- 


10 


18.6 


< 10-10 


-1.074±0.01l3 


0.135±0.015=' 


1.407±0.13l3 


23.7±2.3 


Vela 55(3-12 keV) 


289.50±1.70 


19.3 


1405 


1.20 


< 10- 


10 


42.8 


1 X 10-** 


4.721±0.0033 


0.023±0.0043 


1.164±0.1603 


4.5±0.7 



ITq = MJD 50514.59, Gi and ipi given for the common period of 151.43 d. 
2To = MJD 50514.59. 
^To = MJD 44450.00. 
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